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ABSTRACT 

The precision study of dark matter using weak lensing by large scale structure is 
strongly constrained by the accuracy with which one can measure galaxy shapes. 
Several methods have been devised but none have demonstrated the ability to reach 
the level of precision required by future weak lensing surveys. In this paper we explore 
new avenues to the existing Shapelets approach, combining a priori knowledge of the 
galaxy profile with the power of orthogonal basis function decomposition. This paper 
discusses the new issues raised by this matched filter approach and proposes promising 
alternatives to shape measurement techniques. In particular it appears that the use 
of a matched filter (e.g. Sersic profile) restricted to elliptical radial fitting functions 
resolves several well known Shapelet issues. 
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1 INTRODUCTION 

■ Galaxy shapes provide the unique signature of grav- 
, itational lensing by large scale structure, which has 
' been recognized as a key to the stud y of dark mat- 
ter and dark energy (jMunshi et al.l [20081 ). A limiting fac- 
tor is the accuracy with which one can measure shapes 
dHevmans et al.ll2006l : iMassev et alll2007l : iHoekstra fc JainI 
I2OO8I ). Among the different existing methods, one par- 
ticularly interesting approach is the decomposition of 
galaxy images using basis functions e.g. 'Bernstein & JarvisI 
||2002|) or S hamlets (Rcfrcgicr 2003; Rcfrogicr & Bacon 2003; 
iMassev "fc Refregieii |2005| ) . The strengths of this approach 
rely on the fact that the shape measurement is analytical 
and therefore time efficient as it involves rather small ma- 
trix multiplications. Shapelets decompose an image into a 
linear combination of orthonormal components up to some 
truncation order, and the shape parameters are extracted 
from a least-squares best fit using the recomposed (noise 
free) model. However, the Shapelet type approach suffers 
from a few difficulties: 

(i) The choice of the decomposition truncation order 
is arbitrary. In practice, different lensing groups use radi- 
cally different "optimal" truncation orde rs. Some prefer low 
l|Kuiikenll200^ ). while others prefer high (|Berge et al.ll2008l ). 
although the values for different truncation order could 
be very different. Therefore a constant criterion to mea- 



sure the shape does not appear to be a robust guarantee of 
unbiased shape measurement. 

(ii) "Easy cases" such as large and bright elliptical 
galaxies are poorly fitted. This suggests that a good fit for 
low signal-to-noise galaxies does not necessarily mean that 
the shape has been correctly measured, since it could just 
be buried in the sky noise. This is the overfitting problem. 

(iii) Basis decomposition has too many degrees of free- 
dom for shear measurement, since ideally we are only in- 
terested in two numbers (or six if we include the flexion). 
This is where galaxy morphology and shear measurement 
are clearly two different problems. 

All of those problems have one common origin, namely 
the choice of the zeroth order weight functio n - Gaus- 
sian functions for both Cartesian iR cfrcgicr l2003h and 
Polar Shapelets ( Massey fc Refregien |2005j), as well as 



roiar anapeiets 1 iviassey 6c nerregien izuuoi i. as well as 
iBernstein fc Jarvk 1 20021 ). Unfortunately. Gaussian func- 
tions are poor matches to real galaxy profiles. Ideally, we 
would like the zeroth order to be as close as possible to the 
real profile, and leave to the basis decomposition the task to 
fit departures from this "typical" profile. 

Currently the most promising shape meas urement 
method uses a bayesian mo del fitting approach ( Mill ei~et al-l 
I2OO7I : iKitching et"aLl I2OO8I ). This method does not suffer 
from the same issues as Shapelets, but is limited by the 
strong galaxy profile prior. 

In this paper, we investigate how the change of the 



2 W. Ngan et al. 



weight function affects the basis decomposition method, and 
how it leads naturally to a hybrid method which combines 
Shapelets and fitting techniques. We choose to focus on the 
Sersic profile (hence the term Sersiclets), but our discussion 
can be extended to any profil^E e.g. Moffat profile for ground 
based point spread function (PSF). Section [2] introduces the 
notation and gives a technical description of the new fitting 
functions. Section [3] shows the impact of those fitting func- 
tions on shape fitting and decomposition. Finally, Section |4] 
summarizes our work so far and future possibilities. 

Note that in this paper we choose not to discuss the PSF 
deconvolution. Indeed, the problems we mentioned earlier 
affect equally the measurement of galaxy shapes whether or 
not the galaxies are convolved with a PSF, and Shapelets are 
a popular approach because their Gaussian properties allow 
for very efficient PSF treatment. The PSF deconvolution is- 
sue goes beyond this work because it depends on how the 
PSF is measured and interpolated between stars. Moreover 
the approach developed here could as well be applied to the 
PSF profile measurement separately, and then used later to 
address the deconvolution step through a forward convolu- 
tion m odel fitting method. See for example iKitching et al.i 
ll2008t) . 



2 METHODOLOGY 

2.1 Basis functions in polar coordinates 

In ID, all polynomials Pk{x) of degree k are orthonormal 
with respect to a weight function w{x) if they satisfy 



Pi{x)Pj{x)w(x)dx = 5i 



(1) 



A particular choice of weight function w{x) uniquely de- 
termines the family of polynomials (e.g. for the Cartesian 
Shapelets, a Gaussian weight defines the Ifermite polynomi- 
als). A complete set of polynomials can be useful for decom- 
posing an arbitrary function f{x) as a linear combination of 
basis functions X-n-i^) = Pn{x)[w{x)]^^'^ such that 



f{x) = ^A„Xn(2 



(2) 



In 2D, the basis functions can be represented using po- 
lar coordinates. Following the intuition for Polar Shapelets, 
we separate our basis functions Xmn (r, into the radial 
component Rn(r) and the angular component e""'*. We also 
assume that the weight function u;(r) has no angular depen- 
dence. The 2D basis functions Xmn(f, 0) in polar coordinates 
would be in the form 



X™(r,0) = i?„(r)[t«(r)]i/^e™^ 
The orthonormality is then written as 



(3) 



ra r2 

/ r dr I 
Jo Jo 



Xmn (^5 0)XTn' n' (^1 0) ^mm' ^nn' (4) 



^ While working on the concepts discussed in this paper, we be- 
came aware of similar investigations using exponential and hyper- 
bolic sech functions (Kuijken & van Uitert in prep). 



where * denotes complex conjugate. The orthonor- 
mality requirements for radial and angular parts. 
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RnRniw{r)rdr = and /o^e 
can be satisfied independently. In particular, R„{r) is an 
orthonormal polynomial of degree n with respect to the 
weight function w{r). The integration limit a for the radial 
component will be discussed in Section [2.41 

With an orthogonal and complete set of basis functions, 
an arbitrary image f{r, 0) can then be decomposed into 



(5) 



where the complex basis coefficients Amn satisfy A'^^ — 
A-mn SO that /(r, 0) is wholly real. We will refer rimax as 
"order" in the following sections. 

Readers familiar with Polar Shapelets may notice that 
our radial component here only requires n, and m increases 
in steps of 1. In Polar Shapelets, the radial component re- 
quires both m and n, and m increases in steps of 2. Our 
choice to completely decouple the radial and the angular 
components is of mere convenience, which comes with the 
cost that our set of "basis functions" is no longer complete. 
Consequently, our set of fitting functions cannot decompose 
an arbitrary image. As we shall see in Section [^3] and Section 
12.61 though, the lack of completeness is not a hindrance to 
decomposing images of galaxies fo r weak lensing, as galaxies 
follow Sersic profiles (|Sersid[l968t ) and are not arbitrary in 
general. 



2.2 Weight function 

The basis functions in Shapelets often require high order 
polynomials to describe galaxy shapes accurately because 
galaxies' radial light profiles do not match the weight func- 
tions. Galaxies' lig ht profiles are we ll described by Sersic's 
empirical formula (|Peng et al.ll20o3 ): 



I{r) = I{k) exp[-6A(r/A;)'/^ + bx] 



(6) 



where fc is a scale radius, and A is known as the Sersic index. 
For 0.5 < A < 10, 6a = 2A — 1/3. We use a parameterized 
form of Equation [S] as our weight function: 



w{r) = exp 



-(2A- 



l/X 



(7) 



2.3 Radial component 



The radial component involves a non-trivial computational 
step, as Rn{r) must satisfy the orthonormality require- 
ment described in Section 12.11 We obtain R„{r) by the 
Gram-Schmidt procesfl. Using Dirac notation {Ri\Rj) = 
/ Ri{r)Rj{r)w{r)rdr, the Gram-Schmidt process generates 
each R„{r) by a recurrence relation: 



Ru{r) 



(r7?„_i[ii„_i)\ (i?„_i|i?„_i> 
T' tt: r::: ;- 1 Hn-l rr^ rHn-2 



(i?„-2|-R„-2>' 



where 7?o(r) = 1, and Ri{r) = (r - ji^l^) Ro- After the 
recurrence step, each Rn{r) is individually normalized. 

^ http:/ /mathworld. wolfram. com/Gram- 
SchmidtOrthonormalization.html 
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Figure 1. Orthonormal polynomials of degrees 1, 2, 3, and 4 with 
respect to ■u;(r) = exp [-7.67(r/0.5)^'''']. Orthonormality holds 
in < r < 1 (upper panel) and < r < oo (lower panel). 
The functions in the lower panel are hardly distinguishable, hence 
lack linear independence. This is due to large variation in the 
polynomials coefficients, so the high order details are not visible. 
The functions in the top panel are distinguishably different. 

2.4 Integration limits 

When generating the radial polynomials (Section l2.3|l . a sen- 
sible integration limit must be chosen. In Shapelets one can 
integrate r from to co thanks to the Gaussian function's 
localized profile. For Sersic functions in general, however, 
the profile may not be localized enough to allow for an in- 
finitely large domain. This is generally true for any galaxy 
and stellar profile used as a weight function. 

The problem with using an infinitely large domain is the 
lack of mutual independence among the fitting functions. In 
order to construct a model as a linear combination of fitting 
functions, each fitting function must be distinct so that there 
is no redundancy in their shapes. Figure[T]shows the polyno- 
mials that are generated using a weight function (Equation 
O with A = 4. We find that limiting the orthogonality to 
a finite domain preserves linear independency better than 
extending to an infinitely large domain. 

We conveniently choose < r < 1 for our domain. 
For a square image stamp of IN x 2N pixels, r would be 
normalized to have units of XjN pixels. It also allows us to 
constrain < fc < 1, as the scale radius is always positive, 
and a galaxy should be well captured in a stamp. 

2.5 Completeness 

Although the fitting functions of Sersiclets are indeed mu- 
tually orthonormal, they are not necessarily complete. As 
a result. Equation [5] does not necessarily converge, even at 
high nniax- The basis functions in both Cartesian and Po- 
lar Shapelets are complete, meaning that Equation [5] can 
converge for any arbitrary /(r, 0) as rimax oo. 

The incompleteness of Sersiclets is shown in Figure (2] 
It shows the average difference squared per pixel (Apix^) 
when decomposing a noiseless elliptical object on a 128 x 128 
grid by integrating Equation [5] with Xmn, and exploiting or- 
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Figure 2. Convergence of Equation \5\ represented by the av- 
erage difference squared per pixel (Apix'^) when decomposing a 
high resolution and noiseless elliptical profile. Polar Shapelets are 
clearly more efficient at describing arbitrary profiles. 



thonormality to obtain each Amn- The reconstruction us- 
ing Sersiclets does not improve even as the order increases; 
in fact, the reconstruction becomes slightly worse because 
higher order contributions are as small as the discretization 
error. Clearly, Sersiclets fail to decompose even a simple el- 
liptical object. 

For this reason, we would not attempt the decomposi- 
tion by including higher order fitting functions. In the next 
section we take a different approach - throwing away all cir- 
cularly asymmetric components. The lack of contribution by 
those components is an important rationale of our technique 
to reduce the set of fitting functions. 



2.6 Fitting set reduction 

It is clear from Section [2.51 that the fitting set should be re- 
duced in order to take advantage of the matched filter due to 
the lack of completeness. In the following we reduce our set 
of fitting functions to only the circularly symmetric compo- 
nents (m — 0), and we introduce ellipticities by transforming 
the now perfectly circular model by "scaling" and "rotat- 
ing" , which yield unique v alues of ei and 62. T his is similar 
to the process described in iBernstein fc JarvisI (2002). 

The set of "reduced Sersiclets" has two advantages; 
eliminating the m 7^ components not only cures the over- 
fitting problem, but it also offers a dramatic increase in 
speed as the number of terms in Equation [S] now increases 
like O(nmax) rather than 0{ri^^^). It also provides a direct 
estimate of ei and 62, which are treated as asymmetric scal- 
ing parameters for the fitting function. 

We focus on fitting profiles that are smooth, centrally 
peaked, and elliptical in general. These fitting functions are 
not suitable for studying galaxy morphology, as they can- 
not provide information about a galaxy's detailed structure. 
In weak gravitational lensing studies, however, the details 
in the typical faint images analyzed are dominated by noise 
and should not be fitted. Therefore our reduced fitting set 
offers a natural regularization process which is missing in 
the standard Shapelet approach. Our method is a hybrid of 
Shapelets and fitting techniques, where we do allow some 
decomposition into fitting functions, but those fitting func- 
tions are by construction axisymmetric and therefore pre- 
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de Vaucouleurs exponential generic Sersic 




exponential exponential generic Sersic 



Figure 3. A subset of noisy images for our experiment. Top: Cir- 
cular profiles to test the full fitting set. Bottom: Elliptical profiles 
to test the reduced fitting set. De Vaucouleurs and exponential 
profiles have A = 4 and A = 1, respectively. 

vent isophote mixing (i.e. overfitting) as higher order fitting 
techniques do. 



3 EXPERIMENT 

Our experiment at this stage is not a rigorous test for shape 
measurement, as our test cases (Figure [3| are idealized pro- 
files without PSF convolution. Rather, we are exploring the 
effect of using a variety of weight functions with different 
(fe. A). Our test cases consist of both circular and elliptical 
profiles. We generated two-dimensional reduced x maps of 
k vs order at fixed A values. 

The x?cd maps are then compared against those gen- 
erated using Polar Shapelets, which we will refer to sim- 
ply as "Shapelets" in the following discussion, fn Shapelets' 
case, the r coordinate is also normalized to < r < 1. The 
"scaling-factor" (3 in Shapelets is now comparable to k in 
Sersiclets, which is relative to the size of the image. The 
model fits for both Sersiclets and Shapelets were computed 
using HROTHGAflfl implemented in C. 



3.1 Results — Circular models using the full 
fitting set 

We first test the full fitting set by fitting against circular 
profiles. In Figure |4l we see that convergence of x?od ~ 1 is 
achieved very quickly. This is not surprising as the weight 
functions in the models are indeed realistic. More impor- 
tantly, we find that the fits are insensitive to the choice of 
(fe. A) after the first few orders. This robustness allows us 
to obtain good fits without searching for an optimal (fc,A). 
This result is useful in fitting large collections of objects, 
where families of objects can simply share the same pair of 
(fc. A) without compromise. 

Comparing the x?od maps and image reconstructions 



^ http:/ /hrothgar. sourceforge.net/ 
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Figure 4. Model fits to the de Vaucouleurs (1/A = 0.25) profile 
using the full fitting set. Upper panels: Sersiclets (order 1, 1/A = 
0.3, fc = 0.3) (left) and Shapelets (order 4, A; = /3 = 0.2) (right). 
A^rod ™liiss are 0.9371 and 0.9737, respectively. Lower panels: X^cd 
maps of the same fits as the upper panels. 

between Sersiclets and Shapelets reveals Sersiclets' advan- 
tage. In very non-Gaussian cases such as the de Vau- 
couleurs case, Sersiclets converge at a much lower order than 
Shapelets. In fact, Shapelets require low signal-to-noise ra- 
tios in order to render an illusion of "good fit" . From Figure 
|4l we see that Shapelets' lowest order best fit at order = 4 
already shows signs of noise fitting; the model is not smooth, 
and it shows non-circular isophotes. Sersiclets, though, can 
recover the smooth and circular profile at orders or 1. 

3.2 Results — Elliptical model using the reduced 
fitting set 

For the elliptical profiles, the x^d maps (Figure [S]) are very 
similar to those shown in Figure |4l This means that the m 7^ 
components were indeed not important, and Sersiclets' 
robustness in (fe. A) are preserved. Fitting set reduction as 
described in Section|2]6]has been done to both Sersiclets and 
Shapelets. In Shapelets' case, since m = Q components do 
not exist for odd orders, only even orders were possible. As 
seen in the image reconstruction in FigureO the Shapelet fit 
no longer shows noise fitting thanks to fitting set reduction. 

The residuals of measured ei values corresponding to 
each input value have been plotted in Figure |6l It is clear 
that as the fit improves with higher orders, the scatter in the 
measured ellipticity is reduced. For Shapelets, the solution 
does not show robustness in k as the measured ellipticities 
are more scattered than the Sersiclets' case. 

3.3 Discussion 

Sersiclets' robustness in (fc. A) and its ability to converge 
in relatively low orders comes from the abundance of its 
degrees of freedom in the model. For the full fitting set. 
Equation [S] has (rimax + 1)^ terms in the summation as the 
angular quantum number m increases in steps of 1. Polar 
Shapelets, however, have only (wmax + 1) (wmax + 2) /2 terms 
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Figure 5. Model fits to the elliptical exponential (1/A = 1) profile 
using the reduced fitting set. Upper panels: Sersiclets (order 4, 
1/A = 0.8, k = 0.2) (left) and Shapelets (order 4, k = (3 = 0.2) 
(right). Xrcd are 1.0000 and 1.0634, respectively. Lower 

panels: X?od m^-ps of the same fits as the upper panels. 



Table 1. Input parameters for the elliptical models in Figure |5] 
Each model is 16x 16 pixels in dimension, and has peak S/N 25. 



ei 


62 
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1/A 


0.130 
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0.43 
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0.24 
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0.24 
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0.24 
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0.012 
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0.29 


0.75 
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0.024 


0.29 


0.75 


-0.066 


0.024 


0.29 


0.75 


-0.140 


0.348 


0.28 


0.75 


0.116 


0.357 


0.28 


0.75 


-0.367 


-0.078 


0.28 


0.75 



as m increases in steps of 2 at a given order of n. For the 
reduced fitting set, Sersiclets have rimax + 1 terms, and Po- 
lar Shapelets have only nmax/2 -|- 1 terms. Together with A 
and k, Sersiclets would have about twice as many degrees of 
freedom as Polar Shapelets. To facilitate a fair comparison. 
Figure [6] compares Sersiclets at order 3 against Shapelets at 
order 6 in the second and third panels. This ensures that the 
test uses the same degrees of polynomials available to each 
decomposition technique. 

A drawback of Sersiclets is numerical instability. Her- 
mite or associated Laguerre polynomials in Shapelets can be 
generated very easily with elementary operations, whereas 
the general Sersiclet polynomials require gamma functions. 
In particular, the analytical solution to the integral 

/ r^exp[-&A(r/fc)'/^]dr (i = l,2...) (8) 
Jo 

which is ubiquitous in Section [2.31 can be written as 
X{bxk-^^^)-^-'^ [r(A + jA) - r(A + jA, bxk-^^^)] (9) 
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Figure 6. ei residuals measured using the reduced fitting set. 
Each column of points is an object (Table [TJ, and each data point 
is a diflferent (fc,A), where 0.05 A; ^ 0.95 in steps of 0.1 and 
0.4 ^ 1/A ^ 1.0 in steps of 0.2 (for Sersiclets). The behavior of 
the best fit 62 is very similar to those in ei. At order 3, Sersiclets 
can achieve X^cd ~ ^ measure e with only a small disper- 
sion. In contrast, Shapelets can do so only at a certain range of 
k (or /3) that produces good fits (See Figure [5J. At a given "or- 
der" , Sersiclets have about twice as many degrees of freedom as 
Shapelets. 



where T{z) and r{a,z) are the complete and incomplete 
gamma functions, respectively, and 6a = 2A — 1/3 as before. 
At first glance, it may be tempting to use the analytical so- 
lution because it exists and gamma functions can be readily 
evaluated. Upon closer inspection, though. Equation eval- 
uates the difference between two gamma functions, each on 
the order of (A + jA)!. The difference is then scaled by a 
number raised to the power of —(A 4- jX). One needs to be 
careful when evaluating such an expression if (A -I- j'A)! be- 
comes large, especially when using numerical libraries with 
single precision. Alternatively, the integrand in Equation[S]is 
a smooth function, so the integral can be evaluated directly 
without implementing the analytical solution. 

Sersiclets can be generalized to different weight func- 
tions. In our derivation, although we have chosen the Sersic 
function as our weight function, the process would still be 
the same for any weight functions which are intrinsically el- 
liptical without explicit angular dependence. This allows for 
modeling different types of objects such as the PSF using the 
Moffat profile as the weight function. The integration limit 
of the radial component, as we discussed in Section [2.41 can 
be either finite or infinite depending on whether the weight 
function is sufficiently localized. 
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4 CONCLUSION This paper has been typeset from a TeX/ ]OTb}X file prepared 

• r r,, , , • , • by the author. 

We presented an extension of bhapelets by using an arbi- 
trary weight function in place of the Gaussian function. As 
galaxies' light profiles follow the Sersic profile on average, 
we used the Sersic function as our weight function. This al- 
lowed us to fit cuspy galaxies at lower orders than Shapelets 
could. 

Because the Sersic function lacks analytical properties, 
we used the Gram-Schmidt process to generate the orthonor- 
mal polynomials as radial components for the fitting func- 
tions, where the integrals in the process must be evaluated 
numerically. As the Sersic profile has poor local support, the 
integration limit must be truncated to a finite limit. 

We found that the full set of fitting functions for Ser- 
siclets cannot decompose an arbitrary image even at high 
orders, as the fitting functions do not form a complete set. 
Instead of modeling objects using all fitting functions, we 
reduced the fitting set to only the circularly symmetric com- 
ponents (m = 0). The model was then sheared by ei and 62 
to render elliptical shapes. The reduced set of fitting func- 
tions defines a hybrid method which combines the most in- 
teresting features of the basis decomposition and the fitting 
technique. 

Our experiments so far only focused on idealized im- 
ages simulated using known profiles and noise. Both the full 
and the reduced Sersiclets outperformed Shapelets, as we 
expected. The Shapelet matched filter's true performance 
will be tested in a futur e paper on image simulations such 
as those for GREAT08 (|Bridle et al.ll2008l ). The C code to 
evaluate Sersiclet models is publicly available on request. 
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